function [g,G,Gprime,Gprimeprime,G_generate,weight_function_g,weight_function_b,sup_norm] ...
        = functions_define(method)

%    
% Function for defining the link function, and auxilliary functions
%


if strcmp(method,"atanh")
    g = @(p) atanh(2*p-1); % link function
    G = @(q) (1 + tanh(q)) ./ 2; % inverse of link function p = G(q), q = beta v
    Gprime = @(p) p.^2 + (1-p).^2 ; % derivative of inverse, as a function of p
    Gprimeprime = @(p) 4*p-2; % p-derivative of Gprime
elseif strcmp(method,"logit")
    g = @(p) log(p) - log(1-p); % link function
    G = @(q) exp(q) ./ (1 + exp(q)); % inverse of link function p = G(q), q = beta v
    Gprime = @(p) p.*(1-p) ; % derivative of inverse, as a function of p
    Gprimeprime = @(p) 1-2*p; % p-derivative of Gprime
elseif strcmp(method,"log")
    g = @(p) log(p); % link function
    G = @(q) exp(q); % inverse of link function p = G(q), q = beta' v
    Gprime = @(p) p ; % derivative of inverse, as a function of p
    Gprimeprime = @(p) ones(size(p)); % p-derivative of Gprime
elseif strcmp(method,"linear")
    g = @(p) p; % link function
    G = @(q) q; % inverse of link function p = G(q), q = beta' v
    Gprime = @(p) ones(size(p)); % derivative of inverse, as a function of p
    Gprimeprime = @(p) zeros(size(p)); % p-derivative of Gprime
elseif strcmp(method,"probit")
    g = @(p) sqrt(2) * erfinv(2*p-1); % link function
    G = @(q) (1/2) + (1/2) * erf(q/sqrt(2)); % inverse of link function p = G(q), q = beta' v
    Gprime = @(p) exp(-(erfinv(2*p-1))^2) / sqrt(2*pi); % derivative of inverse, as a function of p
    Gprimeprime = @(p) -2*sqrt(2) * erfinv(2*p-1); % p-derivative of Gprime

end

% Function used to generate data (need not be the same as used in
% Newton-Raphson -- best if maps to interval [0,1]
    G_generate = @(q) (1 + tanh(q)) / 2; % this one based on atahn
    %G_generate = G; % default -- same as model being used
    
% Other functions used in calculations
sup_norm = @(v,w) max(abs(v-w)); % sup norm metric for comparing vectors
%
% The weight functions
%
weight_function_g = @(p,y) 1/(y-1+p); % good weight function (SAMLE)
weight_function_b = @(p,y) (y-p)./(p-p.^2); % bad weight function (standard)

